

clear


%% data and prior
load('../../true_DGP/priorWSigma.mat')
data=readtable('../../../Data/data_main_sample.csv');

data=data(data.year>=1950&data.year<=2000&isfinite(data.D),:);
ccodes=unique(data.ccode);
years=unique(data.year(data.year>1950));
N=length(ccodes);
T=length(years);
K=0;
y=zeros(N,T); D=y; M=y; X=repmat(y,1,1,K);
for n=1:N
    dn=data(data.ccode==ccodes(n),:);
    for t=1:T
        if isempty(setdiff(years(t),dn.year))
            if isfinite(dn.y(years(t)==dn.year))
            	M(n,t)=1; % not missing
                y(n,t)=dn.y(years(t)==dn.year);
                D(n,t)=dn.D(years(t)==dn.year);
            end
        end
    end
end
clear dn t

dd=readtable('../../../Data/distance_capitals.csv');
Z=10^(-6)*table2array(dd(:,2:end));
KZ=size(Z,3);

ep_sd=sqrt(diag(Sigma));

% Sparse-grid integration
[sgi_nodes,sgi_weights]=nwspgr('KPN',2,8);

cnames=cell(N,1);
for n=1:N
    a=unique(data.idacr(data.ccode==ccodes(n)));
    cnames{n}=a{1};
end
ccodes=table(ccodes,cnames);
% World Regions:
% America: 2-165; Europe: 200-255,305,325,375-390,290,310-317,339-370;
% Africa: 404-625,651; 
% Asia & Oceania: 371-373,630-645,652-705,750-771,780-790,710-740,775,800-920
ccodes.region=ccodes.cnames;
for n=1:size(ccodes,1)
    if isempty(setdiff(ccodes.ccodes(n),2:165))
        ccodes.region(n)={'AME'};
    elseif isempty(setdiff(ccodes.ccodes(n),[200:255,290,305,310:317,325,339:370,375:390]))
        ccodes.region(n)={'EUR'};
    elseif isempty(setdiff(ccodes.ccodes(n),[404:625,651]))
        ccodes.region(n)={'AFR'};
    elseif isempty(setdiff(ccodes.ccodes(n),[371:373,630:645,652:705,710:740,750:771,775,780:790,800:920]))
        ccodes.region(n)={'ASO'};
    end
end

clear a cnames data dd n

R=dummyvar(categorical(ccodes.region)); % regions


%% for Knitro

% Jacobian sparsity pattern:
JPattern=[repmat(eye(N),T,1),ones(N*T,2*size(R,2)),ones(N*T,2*N),ones(N*T,N),...
    repmat(eye(N),T,1),zeros(N*T,N),ones(N*T,K),ones(N*T,KZ),eye(N*T)];
JPattern=sparse(JPattern);

% options
options=optimset('Display','off','JacobPattern',JPattern);


%% starting values

% "short" replication
load('thetas_replication.mat','MAP')
pars0=MAP;

% % "long" replication
% load('thetas_replication.mat','pars0')


%% maximum-a-posteriori estimator
[MAP,lp,exit]=knitromatlab(@(x)log_posterior_thetas(x,D,M,X,Z,R,...
    a0,om_a,th0,om_th,b0A,b0D,om_b,s_v,d_v,f0,om_f,s_vs,d_vs),pars0,[],[],...
  [],[],[-Inf(3*N+2*size(R,2),1);zeros(N,1);-Inf(N,1);zeros(N,1);-Inf(K,1);zeros(KZ,1);-Inf(N*T,1)],...
  Inf(length(pars0),1),@(x)eq_con_jac_thetas(x,y,D,M,X,Z,R,sgi_nodes,sgi_weights,ep_sd,Sigma),[],...
  options,'knitro.opt');


%%
clearvars -except MAP exit
save('MAPthetas.mat')



